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ABSTRACT 

Runaway growth ends when the largest protoplanets dominate the dynamics 
of the planetesimal disk; the subsequent self-limiting accretion mode is referred 
to as “oligarchic growth.” Here, we begin by expanding on the existing analytic 
model of the oligarchic growth regime. From this, we derive global estimates 
of the planet formation rate throughout a protoplanetary disk. We hnd that 
a relatively high-mass protoplanetary disk (~ 10 x minimum-mass) is required 
to produce giant planet core-sized bodies (~ 10 Mq) within the lifetime of the 
nebular gas 10 million years). However, an implausibly massive disk is needed 
to produce even an Earth mass at the orbit of Uranus by 10 Myrs. Subsequent 
accretion without the dissipational effect of gas is even slower and less efficient. 

In the limit of non-interacting planetesimals, a reasonable-mass disk is unable 
to produce bodies the size of the Solar System’s two outer giant planets at their 
current locations on any timescale; if collisional damping of planetesimal random 
velocities is sufficiently effective, though, it may be possible for a Uranus/Neptune 
to form in situ in less than the age of the Solar System. We perform numerical 
simulations of oligarchic growth with gas, and hnd that protoplanet growth rates 
agree reasonably well with the analytic model as long as protoplanet masses are 
well below their estimated hnal masses. However, accretion stalls earlier than 
predicted, so that the largest hnal protoplanet masses are smaller than those 
given by the model. Thus the oligarchic growth model, in the form developed 
here, appears to provide an upper limit for the efficiency of giant planet formation. 


Subject headings: Accretion — extrasolar planets — Jovian planets — origin, 
Solar System — planetary formation 


1. Introduction 

The initial growth mode in a disk of accreting planetesimals is runaway growth (eg. 
Wetherill and Stewart 1989, Kokubo and Ida 1996), wherein the mass doubling time for the 
largest bodies is the shortest. However, when these runaway bodies, or protoplanets, become 
sufficiently massive, it is their gravitational scattering (often called viscous stirring) which 
dominates the random velocity evolution of the background planetesimals, rather than the 
interactions among the planetesimals. Since the accretion cross-section of a protoplanet is 
smaller among planetesimals with higher random velocities, protoplanet growth now switches 
to a slower, self-limiting mode, in which the mass ratio of any two protoplanets at adjacent 
locations in the disk approaches unity over time. Ida and Makino (1993) investigated this 
transition analytically and through N-body simulations, and Kokubo and Ida (1998, 2000, 
2002) studied the subsequent accretion mode, giving it the name “oligarchic growth”. In the 
terrestrial region, the hnal accretion phase likely consisted of the merging of oligarchically- 
accreted protoplanets; simulations show that such a process fairly readily produces bodies 
with masses comparable to present-day terrestrial planets (eg. Chambers and Wetherill 
1998). However, an analogous phase in the trans-Saturnian region would have been highly 
inefficient; even sub-Earth mass protoplanets excite each other to high random velocities 
on a timescale short compared to their collision timescale, so that only negligible accretion 
occurs (Levison and Stewart 2001). Thus, it appears that oligarchic growth alone must be 
called upon to account for almost all accretional growth in the outer Solar System. 

In Section 2, we examine the condition for crossover from runaway to oligarchic growth, 
and show that this transition is expected to set in when the largest bodies are still several 
orders of magnitude below an Earth mass. In Section 3, we summarize the previous work 
on oligarchic growth timescale estimates, obtain a protoplanet mass function, then extend 
the model by considering a system in which the planetesimal surface density varies in a 
self-consistent way. We show that with an approximately tenfold increase in surface den¬ 
sity relative to the minimum-mass model, protoplanets of mass ~ 10 M® can form. The 
standard nucleated instability model of gas giant formation, wherein a massive gas envelope 
accumulates onto a solid core during the nebular gas lifetime (^ 10 million years, e.g. Strom, 
Edwards and Skrutskie 1990), is thought to require bodies of this mass (Mizuno et al 1978, 
Pollack et al 1996). Our estimate of the required density enhancement is somewhat higher 
than that of Weidenschilling (1998), who hnds, using a multizone statistical simulation, that 
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4x the minimum mass is insufficient to produce giant planet cores, but that an additional 
“modest increase” is sufficient to make it happen. 

The formation of “ice giant” planets like Uranus and Neptune at stellocentric distances 
of 20 AU cannot be similarly accounted for during this time. In Section 5, we discuss the 
validity of the model. In Section 6, we obtain oligarchic growth rate estimates in the absence 
of gas, to ascertain how much more accretion could have taken place subsequent to the 
removal of the nebular gas. We consider two extremes: that of collisionless planetesimals, 
and that of maximally effective collisional damping of random velocities (though without 
fragmentation). In the former limit, Uranus- and Neptune-mass planets cannot be produced 
at their current locations on any timescale unless the initial protoplanetary disk is implausibly 
massive; in the latter limit, such planets might be formed in a reasonable-mass disk and in 
less than the age of the Solar System. 

In Section 7, we test the semianalytic predictions for the pre-gas dispersal phase of 
oligarchic growth against numerical simulations. We find good agreement as long as proto¬ 
planet masses are well below their theoretical hnal masses, however growth in the simulations 
stalls early, so that the final masses fall short of those predicted by the model. We summarize 
the results and discuss implications in Section 8. 


2. Transition to oligarchic growth 

Ida and Makino (1993) derive the following condition for the dominance of protoplanet- 
planetesimal scattering over planetesimal-planetesimal scattering in determining the random 
velocity evolution of the planetesimal disk: 


2'£mM > 


( 1 ) 


where M and m are the protoplanet mass and the effective planetesimal mass, respectively, 
Sm is the surface mass density of the planetesimal disk, and Sm is the effective surface 
density of a protoplanet in the disk. The latter is given by 
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where a is the semimajor axis of the protoplanet, and Augtir is the width of the annulus 
within which the disk is gravitationally stirred by the protoplanet. For a single protoplanet, 
this width is about 5.2a(e^)^/^, where is the RMS eccentricity of planetesimals in 

the disk. Hereafter, we shorten the notation to e^, and likewise we write 
This width is set by the conservation of the Jacobi energy of planetesimals relative to the 
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protoplanet; details are given in Ida and Makino (1993). A simplification inherent in this 
approach, and one we shall make use of throughout, is the representation of the planetesimals 
as a population of uniform-mass bodies. As discussed by Ida and Makino (1993), the effective 
mass is the RMS mass of the full distribution. This requires that and im depend only 
weakly on m. As we shall see below (Eq. 10), this is indeed the case. At the same time, 
though, one should bear in mind that the planetesimal mass spectrum will span many orders 
of magnitude; as an example, 1 to 100 km radii correspond to a mass range of 10 “^^ to 10 “® 
M 0 . The effect of planetesimal sizes on accretion is considered in Section 4. 


When a disk is stirred by multiple protoplanets, one must account for the overlap of 
stirred annuli. One can do this in an average way by taking 


Aa = Min(5.2ae, Aaproto), 


( 3 ) 


where Aoproto is the characteristic spacing of protoplanet orbits. Kokubo and Ida (1998) find 
that adjacent protoplanets in a swarm of planetesimals keep a typical separation of a fixed 
number b of Hill radii, with b ~ 10. The Hill radius of a body of mass M and semimajor 
axis a orbiting a primary of mass M* is defined as 
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( 4 ) 


where Hm is the reduced Hill radius. Dispersion-dominated random velocities means that 
e 2h, so that lOrn < 5.2ae, and thus we use Aa = lOrn throughout. From Eqs. 1 and 2, 
the protoplanet mass at which oligarchic growth commences is then 
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For the case of an isolated protoplanet (6 ~ 35 — 55) growing in a swarm of m = 10^^ — lO^^g 
(200-400 km radius) planetesimals, Ida and Makino (1993) calculate Mou ~ 10“^ — 10“^ 
M 0 . Considering a population of planetesimals {b ~ 10; see above), and using a perhaps 
more realistic (eg. Lissauer 1987) planetesimal mass of 10^® g (~10 km) one obtains even 
lower crossover masses. Fig. 1 shows Mou throughout the protoplanetary disk for a one 
and ten times minimum-mass nebula (Hayashi 1981; see Eqs. 15 and 17 below) with 10 km 
planetesimals. M^w is a few xl0“^ M 0 or less. Thus, although runaway growth is much more 
rapid than oligarchic growth, it ceases long before protoplanets approaching an Earth mass 
can form. For this reason, the contribution to a planet’s formation timescale from oligarchic 
growth is much more important than that from runaway growth, and we neglect the latter 
in our analysis. 
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3. Oligarchic growth rate estimates 


When planetesimal random velocities are dispersion-dominated rather than shear-dominated, 
the mass accretion rate of an embedded protoplanet is well described by the particle-in-a-box 
approximation (Safronov 1969, Wetherill 1980, Ida and Nakazawa 1989): 


dM 

dt 


j -,2 

~ F——7rR 

h 


M 



( 6 ) 


where h the disk scale height, Rm the protoplanet radius, Uesc the escape velocity from the 
protoplanet’s surface, and Urei the characteristic relative velocity between the protoplanet 
and the planetesimals. F is a factor which compenstates for the underestimation of the 
growth rate which results from using the RMS value of the planetesimal random velocity for 
Urei; its value is ~ 3 (Greenzweig and Lissauer 1992). We apply the following approximations: 
h ~ aim, im — '^esc/'^rei ^ 1) ^^^1 Urei — 6^0^, where hi is the Keplerian frequency. 

These are valid if gravitational focusing and dynamical friction are both effective; see for 
example Kokubo and Ida (1996) for details. The above equation can then be rewritten as 


dM ^ CT.mM^/^ 
dt 

with C = 67r^/^[3/(4pA^)]^^^[G/M*]^/^, where pm is the bulk density of a protoplanet. 


The planetesimals attain an equilibrium RMS eccentricity, e^, when gravitational per¬ 
turbations due to the protoplanets are balanced by dissipation due to gas drag. Following 
Ida and Makino (1993), we obtain by equating the viscous stirring timescale due to a 
protoplanet of mass M, with the eccentricity damping timescale due to gas drag for 

planetesimals of characteristic mass m, The former is given by 
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(Ida 1990, Ida and Makino 1993), where Usm is the surface number density of protoplanets, 
'Em/M. 


The gas drag eccentricity damping timescale is given by 

yem ^ _ ^ ^ I'Q'j 

em (GD/2)7rr^PgasaG 

(Adachi, Hayashi and Nakazawa 1976). Cd is a dimensionless drag coefficient, of order 1 for 
spherical planetesimals of mass 10^® — 10^^ g , is the radius of a planetesimal, and pgas is 
the gas volume density. 
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Setting T-^ ™ and solving for Cm, one obtains 
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where pm is the planetesimal bulk density and the gas density is a power law, pgas oc a ". 
Substituting Eq. 10 into Eq. 7, one gets an estimate for the oligarchic-regime growth 


rate: 


dt 




( 11 ) 


As a check on the validity of using in calculating the accretion rate of proto¬ 

planets, one can compare T®”", the timescale to reach (obtained from Eq. 10 and Eq. 8 
or 9), to the growth timescale, 

T\/r 

( 12 ) 


^ _ M 

= dM/dt 


Using power law gas and solids densities, pgas oc a " and oc a and setting Ud = 1, 
Pm = Pm = 1-5 g/cm^, one obtains 
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For a nebula of ten times the minimum mass or less (see below), this is less than one when 
the protoplanet mass M is a few x 10“^ M® or larger. 

Kokubo and Ida (2000) calculate Tgrow at the heliocentric distances of each of the giant 
planets in the Solar System. But one can also directly solve the above differential equation 
for the protoplanet mass: 
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where we have taken b = 10, pm = Pm = 1-5 g/cm^, Cd=1, and M* = Mq in the 
second line. In this way, one obtains an estimate of the mass of protoplanets throughout 
the disk at any given time. Examples are shown in Fig. 2, which plots protoplanet mass 
versus semimajor axis at various times. The calculation is performed using a surface density 
of solids based on the minimum-mass model of Hayashi (1981): 

= 7.1(a/l AU)-3/2 g/cm^ a < 2.7 AU 
= 30(a/l AU)-/^ g•/cm^ a > 2.7 AU 


where the discontinuity at 2.7 AU is due to the formation of water ice at that heliocentric 
distance (the “snow line”). In reality this boundary may have been significantly less sharply 
dehned. Likewise, the snow line’s location is quite uncertain; other models place it around 6 
AU (Boss 1995) or 1 AU (Sasselov and Lecar 2000), and it was likely not stationary as the 
disk and its temperature profile evolved over time. In our model disk, we spread the snow 
line discontinuity over a radial distance of about 1 AU: 
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For the gas component of the disk, we use the midplane value of the minimum-mass 
density; this is given by 

p““(r) = 1.4 X 10-9(r/l AU)"^^/^ g/cm^ (17) 


The full three-dimensional gas density is of the form 

p(r, z) = po{r) /zoij-f} g/cm^ (18) 

where the disk half-thickness, determined by the balance between the central star’s gravity 
and the gas pressure gradient in the vertical direction, is 

zoir) = 0.0472(r/l AU (19) 


Fig. 2 shows a “front” of protoplanet growth sweeping outward through the disk over 
time. The front is quite steep, due principally to the strong dependence of the protoplanet 
mass on the plantetesimal surface density (Moc S^). For the same reason, the discontinuity 
in the surface density produces a large jump in the protoplanet masses , so that there is in 
effect a second, superimposed growth front which is launched at the snow line. 

Since the surface density is taken as time-invariant in this estimate, the protoplanet 
mass everywhere in the disk increases without bound for increasing time. This is of course 



unphysical; in reality the growth of the protoplanets will be constrained by the total amount 
of solids in the disk. At a given location in the disk, the planetesimal surface density will 
change over time due to depletion of planetesimals by accretion, and due to the systematic 
radial motion of planetesimals caused by nebular gas drag. To obtain an estimate of the 
former effect, we begin by assuming a series of protoplanets spaced by Aa = fern (see Section 
2). One can estimate the planetesimal surface density in an annulus centered on each one: 


0 M 
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where S^ is the original planetesimal surface density at that location. Setting S^ = 0 
also gives us a limiting protoplanet mass, at which growth stops because the surface density 
within Aa drops to zero: 
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This is simply the isolation mass (eg. Lissauer 1987; note that our 6 is a factor of two larger 
than the B in his Eq. 11), modihed as in Kokubo and Ida (1998, 2000) for the case of a 
population of planetesimals, so that b is determined by the orbital spacing of protoplanets 
rather than by the gravitational reach of a single protoplanet. 


Differentiating Eq. 20, one gets the relationship between the surface density and proto¬ 
planet mass rates of change: 


dJ^m _ dM 

accr dt 


( 22 ) 


In the limit of no planetesimal migration, the behaviour of the system is obtained by 
substituting Eq. 20 into Eq. 11: 

— ~AM2/3(S^-5M2/3) (23) 

dt 

with 

Solving this differential equation for the protoplanet mass yields 
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However, to obtain a self-consistent description of the coupled evolution of the proto¬ 
planet masses and the planetesimal disk in the presence of nebular gas, we must also take into 
account the planetesimals’ orbital decay. The rate of change of a planetesimal’s semimajor 
axis under the action of gas drag is given by 


da 

dt 
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^8 2 " 
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(25) 


(Adachi, Hayashi and Nakazawa 1976). Tg^s is dehned in Eq. 9 (note that differs from 
Tgas by a factor of l/cm), a as before gives the exponential a-dependence of the gas density. 


and 


V = 
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is the fractional difference between the gas velocity Vg^s and the local Keplerian velocity uk 
arising from the partial pressure support of the gas disk; Cg is the sound speed, f3 gives the 
exponential a-dependence of the temperature prohle (T oc a“^), and Cs/uk — ZQ{r)/r. Since 
the rate of orbital decay grows with planetesimal eccentricity and inclination, one can expect 
this effect to play an increasingly important role as protoplanets grow larger (Eq. 10). 


Applying continuity, the rate of change of the surface density at a given radius a in the 
planetesimal disk due to planetesimal migration alone is 


dt 
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Using in Eq. 25, we get Vm as a function of a and M{a). We can then 

write the equation of continuity as 
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(28) 


Implicit in this treatment is the assumption that the net orbital decay rate of the planetesi¬ 
mals is not signihcantly affected by the presence of the protoplanets. Tanaka and Ida (1997) 
performed numerical simulations of protoplanets embedded in a swarm of planetesimals, and 
found that a protoplanet shepherds planetesimals outside its orbit, even when mutual per¬ 
turbations among the planetesimals are strong enough to prevent resonant trapping. This 
shepherding effect arises due to the (approximate) conservation of the planetesimals’ Jacobi 
energy relative to the protoplanet. However, they also found that once protoplanets are 
spaced by 15 Hill radii or less, adjacent protoplanets’ contours of constant Jacobi energy 
overlap too much and shepherding no longer occurs. Since Aoproto is closer to 10 rn for a 
population of protoplanets embeddeded in a planetesimal disk (see Section 2), shepherding 
is unlikely to play much of a role. 
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The full oligarchic growth model is thus described by a coupled pair of partial differential 
equations. The hrst is just Eq. 11 with M = M{a,t), while the second is 

the sum of of the two surface density rates of change: 


dE 
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dYjm 

dt 
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dt 
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(29) 


We solve this system numerically, for an initially minimum-mass disk identical to the 
one used to produce Fig. 2. The evolution of M and is plotted in Fig. 3. Protoplanet 
growth now stalls when the planetesimal surface density approaches zero. The results for 
and M in the limit of no planetesimal migration, given by Eqs. 20 and 24, are also shown; the 
latter approaches a limiting value of Miso. The full solution, however, does not quite reach 
Miso because of the depletion of planetesimals by gas drag orbital decay. With or without 
planetesimal migration, by 10 Myrs the largest protoplanet mass is still less than 1 M 0 . 
However, as mentioned in Section 1, the core mass required to form a gas giant planet by 
nucleated instability is likely ~ 10 Mq. In our own Solar System, measurements of the mass, 
radius and gravitational moments of Jupiter and Saturn constrain their present-day solid core 
masses to be 0 - 10 M 0 and 6-17 M 0 respectively (Guillot 1999); it is possible the cores were 
originally larger, but have since become partially mixed with their envelopes (Guillot 2001). 
Assuming a protosolar ice/rock ratio, the core masses of Uranus and Neptune are 85-95% 
of their total masses, that is, 12 - 14 M 0 and 14 - 16 M 0 respectively. This implies that 
the protoplanetary disk out of which our Solar System formed was well above the minimum 
mass, and that a minimum-mass disk is in fact unlikely to produce any giant planets at all. 

Fig. 4 shows the calculation repeated for a disk which has its solids and gas densities 
increased by a factor of ten relative to the minimum-mass disk. This yields a solids density 
of about 25 g/cm^ around 5 AU, which is a litte below the estimate of Lissauer (1987) for 
the density needed to accrete Jupiter’s solid core before the dispersal of the nebular gas. In 
the full numerical solution including planetesimal migration, bodies as large as 10 M 0 form 
just beyond the snow line in less than 1 Myr. Interior to the snow line, the bodies which 
have formed by this time are about an order of magnitude smaller. These results show that 
a strong dependence of protoplanet formation rate on initial solids surface density exists 
even when the depletion of the disk is accounted for, thus lending snpport to the theory that 
the solids density enhancement beyond the snow line facilitates giant planet core formation 
(Morhll 1985). 

In this case, the gas drag-induced orbital decay of planetesimals makes a signihcant 
difference in the later stages of accretion. In the migration-less solntion (dashed curves), 
protoplanet masses get above 60 M 0 ! In the full numerical solution, on the other hand, 
growth already stalls at a little over 10 M 0 due to the rapid orbital decay of planetesimals 
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stirred by protoplanets of this size. 

Also noteworthy is the fact that very little accretion is predicted to take place in the 
trans-Satnrnian region in 10 Myrs. By this time, a lOx minimnm-mass system does not 
prodnce even 1 M 0 bodies beyond abont 15 AU. The protoplanet mass at the location of 
Uranus, ~ 20 AU, is only about a tenth of an Earth masses. The simplest description of 
oligarchic growth, in which the planetesimal surface density is constant (Eq. 14), is a good 
approximation of the early stages of the process, when protoplanets are small compared to 
their isolation mass (Eq. 21). From Eq. 14, the protoplanet mass at a given location and 
time scales as S^. Thus the surface density would have to be higher by at the very least 
a factor of (10/0.02)^/^ = 8, i.e. 80x the minimum mass, in order to produce a body of 
order 10 M 0 at the location of Uranus in 10 Myrs. The minimum-mass protoplanetary disk 
contains a total mass of about 0.02 Mq within 100 AU; 80 times that is 1.6 M 0 , far above 
the typical range of observationally inferred disk masses, ~ 0.01 to 0.2 M 0 (eg. Beckwith 
and Sargent 1996, Chiang et al 2001). Also, such a disk would have a Toomre Q parameter 
less than one beyond about 20 AU, thus being unstable to collapse outside that radius, and 
at best marginally stable throughout most of the region inside (eg. Nelson et al 1998). 

The apparent inability of much of anything to accrete over 10 Myrs in the trans- 
Satnrnian region presents a dehnite obstacle to understanding the formation of our Solar 
System. After all, the 1-3 M 0 gas atmospheres of Uranus and Neptune imply that they are 
composed of protoplanets which, during the lifetime of the nebular gas, grew large enough 
to retain signihcant atmospheres. In fact, the problem of Uranus and Neptune’s formation 
timescale is even more severe, as will be discussed in Section 6. 


4. The effect of planetesimal size 

The efficiency of protoplanet accretion in this model is subject to two competing effects, 
both of which, for a given nebular gas density, depend on the characteristic planetesimal size. 
On the one hand, smaller planetesimals experience stronger damping of random velocities, 
forming a thinner disk and thus increasing the accretion rate. On the other hand, smaller 
planetesimals are also subject to faster orbital decay, which depletes the planetesimal surface 
density at a given location in the disk more rapidly. In other words, a smaller planetesimal 
size means faster accretion, but also an earlier end to accretion. Thus for a hnite available 
time—the lifetime of the nebular gas—there will be an optimal planetesimal radius r™*, 
which produces the largest hnal protoplanet mass. Using, again, the 10 x minimum-mass 
disk, we compute the protoplanet mass after 10 Myrs as a function of planetesimal size at 
different stellocentric distances. The result is shown in Fig. 5. The optimal protoplanet 
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size at 5, 10 and 30 AU is just under 100 km, a few kilometers, and a few tens of meters, 
respectively (one should be a bit cautious about the last value, which is comparable to the 
mean free path of gas molecules at 30 AU for our assumed nebula; this puts us at the 
limit of the validity of the Stokes drag law). Our choice of 10 km planetesimals places us 
logarithmically near the middle of the optimal range for the Jupiter-Saturn region; at the 
same time, even if the planetesimal size at 30 AU had happened to be optimal, the largest 
protoplanet produced there in 10 Myrs would still have been only a few tenths of an Earth 
mass. Another way to look at the above results is in terms of the planetesimal size spectrum 
which will exist in real life (recall that we are approximating the planetesimals as a uniform¬ 
sized population): Accretion is fastest but least efficient from the small end of the spectrum, 
and slowest but most efficient from the large end. 


5. Validity of the estimate 

A number of simplihcations underlie this estimate of protoplanet growth rates. To be¬ 
gin with, interactions among planetesimals are neglected altogether in our analysis. This 
seems reasonable since, by dehnition, scattering by protoplanets dominates the planetesimal 
velocity distribution in the oligarchic regime. Also, Kokubo and Ida (2002) showed that 
the timescale for spreading of the planetesimal disk due to mutual interactions is large com¬ 
pared to the accretion timescale (though the disk does spread signihcantly due to stirring 
by protoplanets, as will be demonstrated by the numerical simulations in Section 7 below). 
However, the protoplanets’ dominance of the planetesimal dynamics also means that the 
relative velocities among planetesimals are large compared to their surface escape velocities. 
Thus, reprocessing of the planetesimal population through physical collisions can potentially 
bring about a signihcant decrease in the characteristic planetesimal size. Inaba and Wetherill 
(2001) performed statistical simulations in the Jupiter-Saturn region; for their adopted frag¬ 
mentation model, they found that a substantial fragmentation tail of small bodies formed, 
resulting in a high loss rate of solids due to gas drag and signihcantly reduced accretion effi¬ 
ciency relative to the case without fragmentation. If fragmentation was indeed very effective, 
it would constitute a signihcant obstacle to the formation of giant protoplanets; it can easily 
be seen from Fig. 5 that if the majority of the planetesimal mass ends up in, say, bodies 
with radii of order 100 m, then the largest protoplanets to form will only be a bit over an 
Earth mass. 

Another ehect of inter-planetesimal collisions will be to dissipate some of the energy in 
planetesimal random motions . We can easily estimate the importance of this ehect relative 
to the damping by gas drag. The timescale for the latter is given by Eq. 9. For the former 
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effect, the lower limit on the timescale is just the time between inter-planetesimal collisions. 
The collision rate can be estimated as n(7rr^)nreh where n is the volume number density of 
planetesimals, n ~ Timlijnh) 'L^/{maem), and Urei CmCiVL. The collision timescale can 
thus be written as 

^^mPm 


T, 


coll 




(30) 


note that it is independent of random velocity. Setting Tcou = and solving for e^, one 
obtains the very simple expression for the eccentricity at which the two timescales are equal: 


2S, 


^coll=gas 


Cy)OjP^ 

as 


(31) 


which is fully determined (within Cd, which we take to be one) by just the gas-to-solids ratio 
plus the prohle of the gas scale height; keeping these hxed, ecoii=gas is unchanged when one 
scales the disk mass up or down. We plot ecoii=gas for a Hayashi-profile disk in Fig. 31. For 
e > ecoii=gas, the gas drag timescale is shorter. Since the collision timescale has to be the 
lower limit to the timescale for random velocity damping by collisions, what we have plotted 
here constitutes the upper limit to the eccentricity at which damping by gas drag comes 
to dominate over damping by collisions. This limit is of order a few times 10“^, thus from 
Eq. 10, inter-planetesimal collisions can only dominate early on, when only small bodies 
have formed and the disk is still dynamically quite cold. It is therefore reasonable to neglect 
collisional damping while gas is present, though one must still remember that collisions may 
reduce the characteristic planetesimal size over time. 


Perhaps most problematic is the issue of migration due to resonant interaction with the 
gas disk. Protoplanets embedded in a gas disk launch density waves at their inner and outer 
Lindblad resonances, and as a result experience a positive and negative torque from the 
inner and outer parts of the disk, respectively (Goldreich and Tremaine 1980). The outer 
torque dominates, resulting in a decay of the protoplanet’s orbit, which has been termed 
Type I migration (Ward 1997). The migration rate increases with protoplanet mass until 
the protoplanet is large enough to open a gap in the gas disk, and then it is locked into the 
slower viscous evolution of the disk (Type II migration). However, since protoplanets likely 
do not form a gap until they reach a mass of ~ 10 - 100 Mq, and since before that, the 
migration timescale at 5 AU can be as short as a few times 10^ years, this poses a problem— 
not just for this approach to estimating growth times, but for our understanding of the 
formation of giant planets in general. It has been proposed that fast inward migration could 
actually speed accretion (Ward 1986); this requires that the protoplanet plough through a 
pristine, dynamically cold planetesimal disk on its way in. However, given the shorter growth 
timescales at smaller stellocentric radii, a migrating protoplanet ought to encounter, instead, 
the dynamically hot remnant of the part of the disk which has already formed protoplanets. 
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Furthermore, simulations have shown that even in the idealized case, accretion efficiency is 
low and the disk has to be enhanced by at least a factor of hve relative to minimum mass 
in order to allow a protoplanet starting at 10 AU to reach gap-opening mass before it falls 
into the star (Tanaka and Ida 1999). 

On the other hand, it should be pointed out that the gravitational interaction of a 
protoplanet with a gas disk is still far from well understood. For one thing, self-consistent 
simulations of multiple non-gap-opening bodies in a gas disk have, to our knowledge, not 
yet been performed; it is not clear what happens when a number of bodies launch density 
waves in close proximity to each other. Conceivably, this could lead to only the inner- and 
outermost members of a population of protoplanets being strongly coupled to the disk. Also, 
Papaloizou and Larwood (2000) hnd that for eccentricities such that the radial excursion is 
equal to the scale height of the gas disk, the direction of orbital migration actually reverses. 
Thus, eccentricity-raising interactions among the protoplanets could in principal counteract 
their orbital decay. However, it is difficult to see how such high protoplanet eccentricities 
(Vr 0.07 at r=5 AU) could arise, since Type I-regime torqnes cause eccentricity to decay 
on an even shorter timescale than semimajor axis. Furthermore, it has been shown (Tanaka, 
Takenchi and Ward 2002) that if the density waves excited by a planet reflect at the outer 
disk edge, the torque asymmetry can be weakened and a nonmigrating steady state may be 
attainable. Finally, it is possible that the gap opening mass is smaller than has been thonght 
thus far (eg. Rafikov 2002). In any case, we defer the issne of Type I migration dnring giant 
planet accretion to future work. For the present, whenever the protoplanet mass somewhere 
in the disk approaches an Earth mass, either in onr analytic estimates or in the subseqnent 
simulations, one should recall that we are dealing with the limit of no disk torqnes. 


6. Oligarchic growth in the absence of gas 

We have established above that collisional damping of random velocities is of little 
importance in the oligarchic growth regime while the gas is present. However, once the 
gas is removed, this may no longer be trne. The issne of post-gas accretion is of particular 
interest in the case of onr Solar System, since from Section 3, the in sitn formation of Uranus 
and Neptnne dnring the gas lifetime appears to be ruled out. Assuming the problem of the 
ice giants’ gas content can be otherwise solved, how well can they grow after the gas is gone? 

Given the limited knowledge which exists abont the behaviour of planetesimals in high¬ 
speed collisions—not a problem amenable to direct laboratory stndy— it is qnite nncertain 
just how effectively these would dissipate energy, and what role fragmentation would play. 
We therefore look at the two extreme cases: a collisionless planetesimal disk, and a disk in 
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which the random velocity damping timescale is eqnal to the collision timescale. 


6.1. The collisionless case 


Withont any damping, an eqnilibrinm planetesimal random velocity no longer exists. It 
thus becomes necessary to simultaneously solve differential equations for the growth rate and 
the evolution of the planetesimal random velocity, v^- Following the approach of Safronov 
(1969), we estimate the latter as 

_ Vm 

dt Trel 



with Trei being the gravitational relaxation time (eg. Binney and Tremaine 1987): 




nM'n'{2GM/vl^Yvm In A’ 


(33) 


where um is the volume number density of perturbers (protoplanets), and A is approximately 
the ratio of the maximum encounter distance (taken to be the disk scale height) to the min¬ 
imum non-collisional encounter distance (taken to be the gravitationally enhanced capture 
radius of the protoplanet); details are given in Stewart and Wetherill (1988). Similar expres¬ 
sions derived for the general velocity evolution of a planetesimal swarm are only valid at low 
random velocities, otherwise the relative velocities between planetesimals are not properly 
accounted for (eg. Stewart and Ida 2000). However, for our simple case of one population of 
bodies on circular orbits stirring another population, taking the relative velocity to be the 
velocity dispersion of the latter should constitute a reasonable (rough) approximation. 

The accretion rate is given by Eq. 6 with Urei = and h = Vm/ Eqs. 6 and 32 
are integrated, using a surface 10 times that of the minimum-mass model. Since the isolation 
mass in the giant planet region of such a massive nebula is about an order of magnitude 
higher than the mass we are trying to attain (~ 10 M®), and since we are interested in the 
upper limit of accretion efficiency, we take to be constant. Protoplanets are taken to 
be initially Mars-mass (0.1 Mq, well above the predicted mass at 20 AU at the time of gas 
dispersal), and the ratio of velocity dispersion to local Keplerian velocity (= e) is initially 
set at 10“^. The accretion rate is set to zero once ^/v‘^ + v'^ is equal to the escape velocity 
from the primary, a/ 2GM*7®, i-®- '^m = vk- Once random velocities have been raised this 
high, the density of planetesimals will have been substantially depleted by ejection, and 
accretion is deemed to have effectively ceased. This is a rather generous estimate of the 
length of the regime over which accretion operates; it can also be argued (eg. Vityazev and 
Perchernikova 1991) that the planetesimal disk ought to already be largely depleted when 
Vm is only ~ 0.3 — 0.4 vk- 
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Fig. 7 shows the results of the post-gas numerical integration at 10® and 10® years. 
Even in 10® years, protoplanets at the locations of Uranus (~ 20 AU) and Neptune (~ 30 
AU) only grow to 4 M 0 and 2 M 0 , respectively. Also plotted is Mstaii, the mass at which 
accretion is throttled by the escape of planetesimals. In the Uranus-Neptune region, it is 
only 4 to 5 M 0 . Thus it appears that in the collisionless limit, a lOx minimum-mass disk 
is unable to oligarchically produce Uranus and Neptune in situ on any timescale. It takes a 
35 X minimum-mass disk to make Mstaw = 10 M 0 at 30 AU. Thirty-hve times the minimum 
mass disk (~ 0.02 Mq) is about 0.7 M 0 . This is several times more than the typical range of 
observed disk masses (~ 0.01 to 0.2 M©; see Section 3). Also, if one supposes that the Solar 
System grew from such a massive protostellar disk, it becomes difficult to explain why it does 
not presently contain more mass—especially in the terrestrial region, where it is difficult for 
material to be ejected. 


6.2. The case of perfect collisional damping 


Next, we look at what happens when planetesimal random velocities are damped on the 
timescale of inter-planetesimal collisions. In this case there exists, once again, an equilibrium 
eccentricity. It is obtained by equating the collision timescale, Eq. 30, to the viscous stirring 
timescale, Eq. 8. The result is 


e 


= 1.7 





(34) 


If we substitute this into the expression for the protoplanet growth rate, Eq. 7, let be 
constant in time and solve, we get 


Mit) 



(35) 


where 

(p,nPMY/hnil^a 

As mentioned in the previous section, the simplifying assumption of a constant planetes¬ 
imal surface density is reasonable as long as protoplanet masses are far below Mnm. In this 
case, since the disk is dissipational, the surface density will also decrease due to planetesimal 
orbital decay. However, since we are interested in the upper limit on protoplanet growth, 
our approximation will suffice. 
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Fig. 8 shows the protoplanet masses computed from the above expression, for the same 
10 X minimum-mass disk as in the previous section. The results for both 10 km and 1 km 
planetesimals are shown. In 10® years, the 10 km case produces neither a Uranus nor a 
Neptune in situ; the 1 km case does exceed the mass of Uranus at 20 AU but only reaches 
a few Earth masses at 30 AU. By a billion years, both cases have exceeded the masses of 
Uranus and Neptune at both 20 and 30 AU. The protoplanet mass is oc as long as 

Sm has not changed much and Mq is small; for example, given that something like a Neptune 
can be grown with = 10 and = 10 km in t ~ 10® years, one can shorten the 
timescale to 10® years by increasing the surface density to a bit under 50 S™™, decreasing 
the planetesimal size to ~ 100 m, or some combination thereof. Thus, it would appear that 
collisional damping does offer a hope for growing large bodies in the trans-Saturnian region. 
However, one must keep in mind that the above growth rate estimates are very much upper 
limits. For comparison, Davis, Farinella and Weidenschilling (1999) report on a statistical 
multizone simulation in the region from 24 to 50 AU, using an approximately 4x minimum- 
mass planetesimal disk. Though they model collisional damping, growth already stalls at 
low masses, and by 4.5 billion years the largest bodies which have accreted are less than 1 
M®. 


6.3. The question of Uranus and Neptune 

The actual growth rate after gas dispersal should, at best, he somewhere between the 
limiting cases of Sections 6.1 and 6.2; even our lower limit may well be overly optimistic. 
Thus, it may not be possible to form Uranus- and Neptune-like planets at all in the trans- 
Saturnian region of a protoplanetary disk, and if it is, the time needed may easily approach 
the age of the Solar System. 

One way to shorten the timescale is to presume a larger disk mass, but as discussed in 
Section 6.1, one can only go so far before coming into conflict with observational results. 
The other possibility is to presume that the planetesimal size in the trans-Saturnian re¬ 
gion was very small, either primordially or as a result of collisional fragmentation, so that 
collisional damping was strong. The primordial size cannot be arbitrarily small, since the 
trans-Saturnian region has to survive millions of years of gas drag without being cleared of 
planetesimals. Fig. 5 gives an idea of the lower limit, since if planetesimals are of the optimal 
size r™* for accretion in the presence of gas (or smaller), this will leave the region largely 
cleared of planetesimals by the time the gas disperses. Fig. 5 shows that is around 20 
m at 30 AU, so during the gas-dominated phase, the characteristic planetesimal size there 
must have been well above that in order to allow the possibility of signihcant post-gas ac- 
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cretion. However, planetesimals could have been collisionally ground down to less than their 
primordial size subseqent to the dispersal of the gas. We defer a more detailed analysis of 
the role of planetesimal fragmentation the post-gas oligarchic growth regime—where it may 
help rather than hinder accretion—to future work. 

An alternative way to account for the existence of Uranus and Neptune in the outer 
Solar System is to lift the requirement that they formed in situ. Zharkov and Kozenko 
(1990) propose that during the final growth phase of a gas giant, it ejects protoplanets 
outward which can serve as the starting point for growing the next giant planet. In this 
way, they suggest that Jupiter triggered the formation of Saturn, which in turn triggered the 
formation of Uranus and Neptune by launching outward protoplanets of a few Earth masses. 
Simulations by Ipatov (1991), performed using a Monte Carlo scheme which neglects distant 
encounters, suggest that this process works if the protoplanets’ eccentricities remain low 
throughout. However, as demonstrated above, even after such a head start, the growth of 
the ice giants could stall well before they reach their present masses. 

The model of Thommes, Duncan and Levison (1999, 2002) assumes that by the time 
Jupiter acquired its massive gas envelope, the Jupiter-Saturn region was able to form, in 
addition to the solid cores of the gas giants, two or more extra bodies of comparable mass. 
They perform simulations which show that the mass increase of a giant protoplanet becom¬ 
ing a gas giant through runaway gas accretion (presumably Jupiter does this hrst) causes 
the remaining protoplanets’ orbits to become unstable. One or more usually undergo close 
encounters with “Jupiter”, and as a result all of them tend to end up on eccentric, mutually 
crossing orbits with aphelia in the trans-Saturnian region. Also included in the simulations 
is a trans-Saturnian planetesimal disk, which serves as a source of dynamical friction for 
the eccentric protoplanets. As a result, the protoplanets’ eccentricities decrease over time, 
decoupling them from Jupiter and from each other on a timescale of a few million years. 
About half of the time, the simulations produce, after 5-10 Myrs, a system which is quite 
similar to our own outer Solar System; Two would-be giant planet cores have ended up on 
nearly circular, low-inclination orbits with semimajor axes similar to those of Uranus and 
Neptune, while the third is near the present orbit of Saturn. The timing of the Saturn core’s 
runaway gas accretion phase is therefore not strongly constrained. Thus, unlike the picture 
proposed by Zharkov and Kozenko, Uranus and Neptune in this model can be transported 
outward to their present locations after they have already completed most or all of their 
growth. Also, scattered protoplanets can be recircularized even after acquiring large eccen¬ 
tricities, by planetesimal disks all the way down to the minimum mass, in contrast to the 
findings of Ipatov (1991). 


This scenario of outward-scattered ice giants appears to £t well with the oligarchic 
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growth model, provided the protoplanetary disk is of sufficient mass to produce the requisite 
numbers and sizes of protoplanets. A disk of 10 times the minimum mass contains 130 
Earth masses of solids between 5 and 10 AU, well above the combined mass of Uranus and 
Neptune, plus the upper limits of the heavy element content of Jupiter and Saturn (Guillot 
1999). As shown in Section 3, bodies of mass ~ 10 M 0 are predicted to form in less than a 
million years. A spacing of ~ 10 th allows about four such bodies to fit in the Jupiter-Saturn 
region. However, the question of how readily this many objects can be produced in the time 
available—while both gas and planetesimals are still present—must ultimately be addressed 
with numerical simulations. 


7. Numerical simulations 

A simple semi-analytic estimate for protoplanet mass as a function of time throughout a 
protoplanetary disk is a potentially powerful tool, since it offers the possibility of character¬ 
izing accretional evolution over time and distance scales which are as yet beyond the reach 
of numerical simulation. Nevertheless, to assess the validity of such an estimate, compar¬ 
isons to simulations must be made. The limits of computing capacity restrict the domains 
of full N-body simulations to relatively narrow annuli within a protoplanetary disk. In the 
simulations presented below, we are able to simulate larger regions of the disk by making 
use of simplihcations which speed computation but preserve enough of the relevant physics 
to make the results meaningful. 


7.1. Method 

The simulations are performed with a variant of SyMBA, a symplectic integrator which 
makes use of an adaptive timestep to resolve close encounters among bodies (Duncan, Levison 
and Lee 1998). This version also models the aerodynamic drag force on planetesimals due 
to a gas disk. The gas disk is modeled with a three-dimensional density prohle of the form 
given by Eqs. 18 and 19, and the parameter rj is calculated from Eq. 26. 

In the vertical direction, the reference Keplerian velocity for the gas disk changes, since 
it is the horizontal component of the Solar gravity which provides the central force: 

v^ia^z) =- —(a/V -f 2 ;^)^ (36) 

a 

Stokes drag is applied to planetesimals: 


V = -KVre\Vre\ 


( 37 ) 
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where Vrei is the velocity of a planetesimal relative to the gas disk. The drag parameter K is 

K = (38) 

^PmXna 

and we again adopt Cjy = 1. 

A number of simplihcations are made to render feasible the task of simulating proto¬ 
planet growth over a radial range of up to tens of AU for millions of years. First, to prevent 
the computational expense from being prohibitive, the planetesimal disks are built up of 
bodies much larger than a realistic characteristic planetesimal (~ 1 to 100 km, correspond¬ 
ing to 10“^^ to 10“® M 0 ). We adopt planetesimal masses of 0.01 to 0.05 M 0 . However, for 
the purpose of calculating the gas drag, the actual planetesimal size is not used. Using the 
approach of Beauge, Aarseth and Ferraz-Mello (1994) we instead assume a more physically 
realistic size of 1 or 10 km (~ 10“^^ or 10“® M 0 ) and apply the drag accordingly. Thus, 
each small body in the simulation can be regarded as a “super-planetesimal” representing 
the averaged orbits of a large number of real planetesimals. 

With such large planetesimals, one must ensure that the protoplanets are sufficiently 
larger, otherwise they will not experience effective dynamical friction. From test runs, we de¬ 
termined that an order of magnitude difference between the populations is enough to keep the 
protoplanets’ eccentricities and inclinations reasonably low. Also, using super-planetesimals 
increases the planetesimal accretion rate of a protoplanet, since the physical interaction ra¬ 
dius is the sum of both bodies’ radii, and since in the realistic case the planetesimal radius is 
negligibly small compared to that of the protoplanet. A factor of ten mass difference gives a 
factor of about two difference in radii, and thus a growth rate oc + ?"m)^ which is initially 
too high by a factor of roughly two. The fractional error decreases as the protoplanets grow 
and widen the protoplanet-planetesimal size gap. 

The large starting mass of the protoplanets constitutes an unrealistic initial condition; 
for example, from the semianalytic estimate, protoplanet masses at 10 AU need several 
million years to reach 0.1 M 0 if they are initially 0.1 M 0 (Fig. 4). Thus the simulations 
give the protoplanets a significant head start. Again, this leads to a smaller fractional error 
in protoplanet growth rates for larger protoplanets, and since the timescale to reach the hnal 
mass is dominated by the time spent in the later phases, the oversized protoplanets are a 
reasonable initial condition. 

Finally, just as in the semi-analytic estimate, planetesimal-planetesimal interactions— 
both gravitational and collisional—are neglected. The planetesimals are treated as a non¬ 
self-interacting population by the integrator, though each one fully interacts with all proto¬ 
planets, thus making these “N-l-N'-body” simulations in which the computation time scales 
quadratically with N (the number of protoplanets) but only linearly with N' (the number 
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of planetesimals). Another benefit of neglecting gravitational interactions among planetesi- 
mals is that it prevents self-stirring of the planetesimal population which, given their large 
masses and lack of softening—at present SyMBA does not support softened potentials due 
to the way close encounters are handled—would result in unrealistically high eccentricity 
and inclination growth rates. (Partially) inelastic collisions between planetesimals, on the 
other hand, would act to reduce their random velocities, but as shown in Section 5, this 
effect is neglible compared to damping by gas drag. In the simulations, collisions between 
protoplanets and planetesimals, or between two protoplanets, are treated as perfectly in¬ 
elastic; that is, the two participating bodies are always merged. Work done on the role of 
fragmentation in late-stage planetary formation (Alexander and Agnor 1998) suggests that 
this is a reasonable assumption. 


7.2. Simulation results 

7.2.1. Run A 

The first simulation is performed in the vicinity of the snow line (2.7 AU) in the Hayashi 
disk. The gas and solids surface densities are increased everywhere by a factor of hve relative 
to the minimum mass model. Equal-mass, 0.01 M 0 planetesimals are initially distributed 
between 1.5 and 5 AU. Planetesimals initially have a Rayleigh distribution in eccentricities 
and inclinations, with RMS values of 0.01 and 0.005 (= 0.29°) respectively, somewhat lower 
than what is given by Eq. 10 (e^ — 0.03 ~ 2im). However, since the timescale to reach 
once the simulation commences is short compared to that for accretion (see Eq. 13), 
unrealistically low initial planetesimal random velocities have little effect on the outcome 
of the simulation. The protoplanets are given equal masses of 0.1 M 0 , and are distributed 
over the same range in semimajor axis, with succesive bodies spaced about 10 rn apart. 
Planetesimal densities are set at 3 g/cm^. The density of ice-enhanced material beyond 
the snow line would have been lower—perhaps only half of this—but since growth rates are 
not a strong function of the body densities, a single density is used for the planetesimals. 
Protoplanets outside the snow line are given an initial density of 1.5 g/cm^, but because 
SyMBA averages densities when bodies merge, their densities also approach 3 g/cm^ as 
they accrete planetesimals. Gas drag commensurate with a size of 10 km is applied to the 
planetesimals. The fractional difference 77 between Keplerian velocity and the gas orbital 
velocity is of order 10 “^. 

The total simulation time is 1 Myr. Snapshots of the protoplanet masses and semimajor 
axes at .05, .1, .5 and 1 Myrs are shown in Fig. 9. Also plotted are the estimates of 
protoplanet mass, obtained from solving Eqs. 11 and 29 with an initial mass of 0.1 M 0 , 



and with the physical protoplanet cross-section increased to take into acconnt the large 
planetesimal radii {tm + rm)- Protoplanet growth stalls earlier than predicted, with 

little change in protoplanet massess between 0.5 and 1 Myrs; there is primarily jnst an 
overall spreading of protoplanet orbits over that time. Since the planetesimals spread with 
the protoplanets, the expansion of the protoplanet system beyond the original radial extend 
of the planetesimal disk acts to lower the planetesimal snrface density, and thns slow growth, 
relative to the analytic description. The largest protoplanet at 1 Myr has a mass of 2 M 0 ; 
the largest model-predicted mass at this time is aronnd 6 M 0 . The assnmption of a 10 th 
orbital spacing is ronghly borne ont by the largest protoplanets. 

A jnmp in protoplanet mass does seem to exist at the snow line for a while, particularly 
at 0.1 Myrs. At later times this jump becomes washed out. But this is to be expected, 
since both protoplanets and planetesimals diffuse in semimajor axis as they gravitationally 
interact with each other, and since the simulation does not “enforce” the snow line once it 
starts running. In other words, sublimation/freezing of ice/water crossing the boundary is 
not modeled. 

Fig. 10 shows the state of the simulation at 1 Myr in more detail, including planetesimal 
eccentricities and inclinations. The eccentricities and inclinations are comparable to the 
predicted values beyond 3 AU. However, inside 3 AU, only two protoplanets remain and 
these have opened a shared gap, thus halting accretion. As a result, the planetesimal random 
velocities here are signihcantly lower than predicted. The effect of migration of planetesimals 
due to gas drag can be clearly seen here. Over time, more and more are deposited interior 
to (most of) the protoplanets, forming a broadening, dense, dynamically cold ring. This 
happens because less eccentric planetesimal orbits decay less rapidly (see Eq. 25), so that 
once planetesimals are out of range of strong gravitational stirring by the protoplanets, their 
radial motion slows and they pile up. Orbital repulsion between the protoplanets and these 
massive rings—at 10® years, there are about 14 M 0 in planetesimals between 2.5 and 3 
AU, and about 16 M 0 interior to 2 AU— further promotes the segregation between the 
planetesimal and protoplanet populations. 

The formation of dynamically cold planetesimal rings like those seen in this and sub¬ 
sequent runs is a simulation artifact, unlikely to occur in a real system. In actuality these 
planetesimals would have been further scattered by protoplanets already formed at smaller 
stellocentric distances. Alternatively, in cases where a high density of planetesimals really 
did accumulate in a region largely devoid of gravitational stirring by protoplanets, they 
would simply revert temporarily to runaway growth and spawn new protoplanets in their 
midst. Therefore, the end state of Run A is almost certainly a poor representation of reality 
interior to 3 AU. In particular, the dense planetesimal ring between 2.5 and 3 AU would 
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have produced much more protoplanet growth, and thus the analytic estimate of protoplanet 
mass would have been better reproduced in that region than it is in the simulation. At the 
same time, as long as a region of the disk is, for whatever reason, predominantly stirred by 
only one or a few protoplanets, the opening of a gap like that in Run A (and subsequent 
runs) is a likely outcome (Rahkov 2001). 


1.2.2. Run B 

The next simulation we present uses a planetesimal disk with an initial surface density 

= 250(a/l AU)“^g/cm^ = 10(a/5 AU)“^g/cm^, (39) 

extending from 5 to 15 AU. This is suggested by Pollack et al (1996) as the optimum surface 
density for the accretion of the giant planets. It is 3.7 times the minimum-mass surface 
density at 5 AU. The corresponding gas nebula is therefore given a midplane volume density 
of 3.7 times its value at 5 AU and a density profile oc under the assumption 

that the gas scale height also has the prohle of Eq. 19. The gas volume density is thus 

= 1.15 X 10-®(a/l AU)-^^/Vcm^ (40) 

The planetesimal disk consists of equal mass bodies of 0.02 Mq. In this disk, protoplanets of 
initial mass 0.2 Mq are placed at intervals of approximately ten Hill radii. For the application 
of gas drag, the planetesimals are, again, assumed to have a size of 10 km. 

The length of the run is 10 Myrs. Snapshots of the protoplanet masses and semimajor 
axes at 0.5, 1, 5 and 10 Myrs are shown in Fig. 11. Fig. 12 gives the state of the run at 10 
Myrs in more detail, showing also the eccentricities and inclinations of the protoplanets and 
planetesimals. The largest protoplanet which has formed by this time has a mass of 2 M0. 
The location of the growth front matches that predicted by the model quite well throughout 
the run. Agreement with the theoretical protoplanet masses beyond 5 AU is likewise quite 
good initially, however the masses reached after the growth front has largely swept past, 
inside 6-10 AU at 1 - 10 Myrs, are smaller than the predicted hnal masses by a factor 
of several. This appears to be at least partly an edge effect: Under the action of stirring 
by protoplanets, the planetesimal disk spreads—this is superimposed on the net inward 
migration of planetesimals—which further lowers the disk surface density, particularly near 
the edges. Since protoplanets grow fastest and orbital times are shortest at the inner edge, 
the surface density there is affected most strongly. Fig. 13 shows that the surface density 
near the inner edge of the protoplanet population, around 5 AU, is indeed lower than the 
analytic estimate after about a million years. As time goes on, this underdense region moves 
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outward; by 10 Myrs, it reaches all the way out to about 11 AU. In addition, Fig. 12 shows 
that, similar to the case of Run A, the innermost protoplanet at 10 Myrs has cut short its 
growth by forming a gap in the planetesimal disk. 

Fig. 14 shows the evolution of protoplanet semimajor axes in detail, including their 
merger history. A total of hve mergers among protoplanets take place, reducing their number 
from the original 16 to 11; no protoplanets are lost to ejection, or otherwise leave simulation 
domain. In this run, as in the others, mergers among protoplanets only play a secondary role 
in the accretion process; sweep-up of planetesimals is the primary growth mode, therefore 
the oligarchic growth model provides a reasonable £t. Also apparent in Fig. 14 is that 
significant reordering of protoplanet orbits takes place over the course of the run. 


1.2.3. Run C 

We repeat the above simulation, but with the drag force on the planetesimals increased 
to what 1 km objects would experience. Snapshots of the run are shown in Fig. 15, and 
the endstate, at 10 Myrs, is shown in Fig. 16. In keeping with what the model predicts, 
the growth front moves outward more rapidly than in Run B. The largest protoplanet at 10 
Myrs is just under 2 M®, very similar to Run B, although the predicted largest protoplanet 
mass is only about half that of Run B. The theoretically predicted discrepancy comes about 
because smaller planetesimals ought to bring about faster growth, at the expense of a smaller 
hnal mass, as discussed in Section 4. However, the analytic protoplanet masses in Run B are 
only larger than those of Run C inside about 8 AU at 10 Myrs. In both runs, the growth of 
protoplanets in the inner region is stunted by simulation edge effects, as described in Section 
7.2.2 above. 


7 . 2 . 4 . Run D 

Lastly, we show the result of a simulation performed with a more massive initial proto¬ 
planetary disk. The solids surface density is chosen as 

=20(a/5AU)-3/Vcm^ (41) 

which is about 7.5 times that of the minimum-mass nebula. The gas surface density is 
increased above minimum mass by the same factor, and the half-thickness is assumed to be, 
again, given by Eq. 19. The total mass in planetesimals is about 2.7 times as large as in 
Run B, so to keep the number of bodies reasonably low, more massive bodies are used in 
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Run D. The planetetesimal mass is changed to 0.05 M 0 , and the protoplanet initial mass to 
0.5 M®. 

The length of this run is 5 Myrs; snapshots at 0.1, 0.5, 1 and 5 Myrs are shown in Fig. 
17, and the endstate is shown in more detail in Fig. 18. The run produces a largest body of 
5 M® just inside 9 AU by 5 Myrs, and a total of four bodies that are more massive than 3 
M®. The initial growth rate is again faster than predicted, but the masses produced once the 
wave of growth has swept past, inside 9 AU at 5 Myrs, are lower than the theoretical value 
by a factor of several. Fig. 18 shows only a small number of planetesimals remaining among 
the protoplanets interior to the growth front, conhrming that accretion there has essentially 
concluded. 


8. Conclusions 

Runaway growth allows very short formation times, but only in the early stages of 
planetesimal accretion; there is a transition to the self-limiting oligarchic growth mode when 
the largest bodies are still orders of magnitude below an Earth mass. The timescale of 
oligarchic growth thus dominates over that of runaway growth, and we use the former alone 
to obtain a global picture of planet formation throughout a protoplanetary disk. In the 
terrestrial region, accretion efficiency is high, and oligarchic growth alone is not required to 
form hnal bodies of terrestrial-planet masses. Simulations show that a late stage of impacts 
among protoplanets of up to perhaps Mars mass is readily able to produce a hnal system with 
planet masses and spacings similar to the present-day inner Solar System (eg. Chambers 
and Wetherill 1998). A signihcant planetesimal population is not necessary to facilitate 
late-stage accretion, though it may be needed to reproduce (via dynamical friction) the low 
eccentricities of the terrestrial planets. 

In the giant planet region, however, lower accretion efficiency and the necessity of form¬ 
ing large bodies before the removal of the nebular gas (after ^ 10 Myrs) mean that direct 
oligarchic growth likely has to be relied on to account for most of the accretion. In our semi- 
analytic model, the outward-sweeping front of oligarchic growth produces 10 M® bodies in 
a protoplanetary disk with ten times the mass of the Hayashi (1981) minimum-mass model 
in less than a million years. These hrst appear near the snow line, where the surface density 
enhancement gives accretion a head start and increases the maximum attainable masses. 
Thus, the snow line may indeed play a role in triggering the formation of giant planets, as 
has been previously suggested (Morhll 1985). 

The numerical simulations performed here agree reasonably well with the model insofar 



as the timing of the initial oligarchic growth front is concerned. However, accretion stalls 
earlier than predicted, and the largest hnal masses produced typically fall short of the the¬ 
oretical values by factors of several. This can be attributed in part to a simulation edge 
effect, namely, the spreading of the simulated planetesimal disk beyond its initial radial ex¬ 
tent, which causes the disk surface density to decrease faster than the model (which does 
not incorporate this effect) predicts; the effect is strongest at the inner edge of the disk. 
Another edge effect is that the innermost few protoplanets act more like isolated bodies in a 
planetesimal disk, because there are fewer neighbouring stirred regions which overlap theirs. 
Consequently, they tend to open a gap in the planetesimals, thus putting a premature end to 
their growth (Rahkov 2001). Gap formation about the innermost protoplanets is evident in 
Figs. 10, 12 and 16. In order to make the inner edge less artihcial, future simulations ought 
to include additional larger protoplanets interior to the planetesimal disk, to play the role 
of the “endproducts” in the region where oligarchic growth has already reached completion. 
However, given that protoplanets interior to the snow line are expected to reach signihcantly 
smaller hnal masses, edge effects like those in our simulations may not be completely unre¬ 
alistic. In any case, it appears that the semi-analytic estimate we have developed here for 
the oligarchic growth endproduct masses is an upper limit. 

Overall, then, it is still somewhat of a challenge to understand how (potential) gas giant 
cores can form in a protoplanetary disk. The formation of an extended gas atmosphere on 
larger protoplanets, not modeled here, will certainly help in sweeping up planetesimals. It 
may simply be that the protoplanetary disks which produce giant planet-bearing systems 
really are quite massive, perhaps in excess of ten times the minimum mass. However, one can 
only push the disk mass so far before one comes up against observational limits. Stevenson 
and Lunine (1988) develop a model in which diffusive redistribution of water vapor from the 
inner part of the nebula leads to a large local density enhancement near the snow line, so 
that the jump in solids density there has an additional spike superimposed on it. Another 
possible scenario is that runaway gas accretion already takes place at smaller core masses, 
perhaps as little as 1 M 0 . This would be facilitated by low grain opacites in the accreting 
atmosphere (eg. Inaba and Wetherill 2001). Yet another possibility is that the window of 
opportunity for gas giant formation is larger than normally assumed. Kokubo and Ida (2002) 
adopt 10® rather than 10^ years as the gas lifetime, pointing out that the timescale for the 
complete dispersal of gas may approach this larger value (Thi et al 2001). However, as the 
gas density drops, so will the planetesimal accretion rate. Also, it is unclear how much the 
gas nebula can be depleted before core accretion stops being viable. 

In 10 Myrs, the (optimistic) predicted growth front only produces 10 M 0 bodies interior 
to about 10 AU, even in very massive protoplanetary disks. Thus in the case of our Solar 
System, although the solid cores of Jupiter and Saturn can be more or less accounted for, it 



seems that Uranus and Neptune are out of luck as far as in situ formation during the lifetime 
of (most of) the gas is concerned. Accounting for the ice giants’ gas content is a problem, 
since the protoplanets which are predicted to have formed in the trans-Saturnian region by 
the time the gas disperses are too small to have acquired appreciable atmospheres. The 
most fundamental difficulty, however, is that of the post-gas growth timescales. If random 
velocities in the planetesimal disk are not damped at all, then the formation of Uranus- and 
Neptune-mass planets is not possible on any timescale. If collisions among planetesimals 
are an effective dissipational mechanism, then such planets could conceivably form, though 
it could well take of order a billion years, during which time the planetesimal disk has to 
maintain a sufficiently high optical depth to keep the collisions going. 

A model in which the ice giants shared the same birthplace as Jupiter and Saturn, only 
to be scattered outward when Jupiter accreted its massive gas envelope, provides an alter¬ 
native (Thommes, Duncan and Levison 1999, 2002). Of course, this scenario is predicated 
on the ability to form multiple ^ 10 Mq bodies in the Jupiter-Saturn zone, notwithstanding 
the above difficulties, while at the same time leaving some of them without massive gas 
envelopes at the time when the nebular gas disappears. The simulations performed here 
provide some support for this scenario, in the sense that they show several comparable-mass 
largest protoplanets occupying the gas giant region after several million years of evolution. 
Run D, with the largest initial disk mass, looks the most promising. However, the proto¬ 
planets produced still fall signihcantly short of Uranus and Neptune’s core masses. Run D’s 
protoplanet masses at 10 Myrs would need to be scaled up by a factor of about three; then 
the largest protoplanets would comprise a system not dissimilar to the initial conditions used 
in Thommes, Duncan and Levison (1999, 2002) 

Orbital evolution through disk tides (eg. Ward 1997) is not considered in this analysis. 
Fast Type I inward migration, predicted to be strongest for bodies with masses of order 10 
- 100 M 0 , constitutes a major potential problem for any model of giant planet formation, 
except ones which invoke direct, unnucleated collapse from a gas disk instability (eg. Boss 
1998). However, at present it is still quite uncertain just how serious this problem is. For 
instance, it has been shown that the nature of the migration depends sensitively on the as¬ 
sumed disk boundary conditions (Tanaka, Takeuchi and Ward 2002). Also, Type I migration 
may be halted at a smaller mass than is commonly assumed (eg. Rahkov 2002). 

In the limit of no Type I migration, the investigation presented here suggests the fol¬ 
lowing overall picture of the post-runaway planet formation process: During the first ten 
million years or less, while the nebular gas is still present, a relatively rapid front of oli¬ 
garchic growth sweeps outward through the disk. By the time the gas disperses, however, 
this wave has only reached what corresponds to the Jupiter-Saturn region; very little accre- 



tion takes place during this time in the trans-Saturnian region, thus the planetary system is 
initially quite compact, with a Kuiper belt starting at ~ 10 AU. Bodies large enough to be 
potential gas giant cores form only in an annulus encompassing roughly the outer half of the 
(proto-) planetary system’s radial extent. Once the gas disappears, the outward-expanding 
wave of growth becomes far slower and the subsequent formation of giant planets (of ne¬ 
cessity gas-poor ice giants like Uranus and Neptune) becomes difficult, perhaps impossible. 
However, if gas giants are able to form before the removal of the gas, a likely by-product will 
be the scattering of any remaining giant protoplanets that missed out on acquiring a massive 
gas envelope. This could ultimately produce an outer planetary system with widely-spaced, 
circular giant planet orbits like those of the present-day Solar System, while requiring little 
or no accretion to take place after the gas is gone. 
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Fig. 1.— Threshold protoplanet mass for crossover from runaway to oligarchic growth in 
a minimum-mass (lower curve) and 10 x minimum-mass (upper curve) nebula, as obtained 
from Eq. 5. A uniform planetesimal size of 10 km is assumed. 
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Fig. 2.— The protoplanet mass function for a minimum-mass nebula at 0.5, 1, 5 and 10 
Myrs. The curves are computed using Eq. 14, which assumes a planetesimal surface density 
that is constant over time, so that the protoplanet mass (unphysically) increases without 
bound. Mq is taken to be zero; as long as the initial mass is well below the mass range of 
interest, it has negligible effect on the solution. 
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Fig. 3.— Protoplanet mass (thick curves) and planetesimal surface density (thin curves) ver¬ 
sus heliocentric distance for a minimum-mass gas nebula (Eq. 17) and an initially minimum- 
mass planetesimal disk with a smoothed-out snow line discontinuity (Eq. 16) after 0.5 to 10 
Myrs of evolution. The solution obtained from numerically integrating the coupled system 
of Eqs. 11 and 29 is shown as solid curves, while the analytic result when migration of 
planetesimals is neglected, given by Eqs. 20 and 24, is shown as dotted curves. The mi¬ 
grationless solution approaches the isolation mass Mum, shown as a thick dashed line, over 
time. A uniform planetesimal radius of 10 km is assumed. 
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Fig. 4.— Protoplanet mass (thick curves) and planetesimal surface density (thin curves) 
versus heliocentric distance for a nebula with 10 times the solids and gas densities of the 
minimum-mass nebula. All other parameters are the same as for the minimum-mass case 
shown in Fig. 3. The dashed curves show the migration-less solution computed from Eqs. 20 
and 24, while the solid curves show the solution including planetesimal migration, obtained 
from numerically integrating Eqs. 11 and 29. 
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Fig. 5.— The largest protoplanet mass after 10 Myrs of oligarchic growth as a function of 
planetesimal radius at 5, 10 and 30 AU, in a 10 x minimum-mass disk. As can be seen, the 
optimal planetesimal radius r"'* is about 80 km, 2 km and 20 m, respectively. 
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Fig. 6.— For a Hayashi-profile disk, the eccentricity ecoii=gas at which the timescale for 
random velocity damping by gas drag comes to be equal to the planetesimal-planetesimal 
collision time. At higher eccentricities, the gas timescale is shorter. 





a(AU) 

Fig. 7.— Protoplanet mass versus semimajor axis after 10® and 10® years of post-gas dispersal 
accretion, in a disk having ten times the solids surface density of the minimum-mass nebula. 
The lower of the dashed curves shows Mstaii, the protoplanet mass at which signihcant escape 
of planetesimals likely terminates accretion. The upper dashed curve shows Mgtaii for a 35x 
minimum-mass disk. 
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Fig. 8.— The protoplanet mass at 10® and 10® years, in the limit of maximally effective 
collisional damping of planetesimal random velocities. The dashed curve shows the limiting 
mass (Eq. 21), which is simply the mass at which a protoplanet has accreted all solids in an 
annulus of width 10 th. The top panel shows the results for 10 km planetesimals, while the 
bottom panel shows the 1 km case; M oc l/r^ so in the latter case the masses are simply 
scaled up by a factor of ten. 
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Fig. 9.— Snapshots of the protoplanet masses and semimajor axes in Run A at different 
times. The superimposed curve shows the mass predicted by the oligarchic growth model at 
each time, and the horizontal “error bars” on each protoplanet have a length of 10 Hill radii. 
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Fig. 10.— A more detailed snapshot of the state of Run A at its endpoint of 1 Myr, showing 
also the planetesimals. Eccentricities (top), inclinations (middle), and protoplanet masses 
are plotted versus semimajor axis, with superimposed curves showing their model-predicted 
values. 
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Fig. 11.— Snapshots of the protoplanet masses and semimajor axes in Run B at differ¬ 
ent times, together with the model-predicted mass function. The horizontal bars on the 
protoplanets are of length 10 th. 
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Fig. 12.— A more detailed snapshot of the state of Run B at its endpoint of 10 Myrs, 
showing eccentricities (top), inclinations (middle) and protoplanet masses (bottom) together 
with their predicted values. 
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Fig. 13.— Evolution of the planetesimal disk surface density in Run B (histograms), reveal¬ 
ing simulation edge artifacts. The superimposed curves show the model-predicted surface 
densities. By about a million years, the planetesimal surface density in the simulation falls 
well short of the theoretical value near the inner edge of the protoplanet population at 5 
AU; by 10 Myrs, the underdense region has moved all the way out to ~ 11 AU. At the same 
time, inward-migrating planetesimals pile up interior to the protoplanets, leading to locally 
very high surface densities. 
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Fig. 14.— Evolution of protoplanet semimajor axes in Run B. The plotted circles are pro¬ 
portional to the protoplanets’ physical sizes. The connecting lines trace the history of each 
protoplanet. Going from bottom up, wherever multiple lines go in and one line comes ont 
of a protoplanet, mergers have occurred in the intervening time; there are hve mergers over 
the course of this run (note that there are no mergers at 10® yrs, only two protoplanets with 
nearly the same semimajor axis, abont 6.5 AU) 
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Fig. 15.— Snapshots of the protoplanet masses and semimajor axes in Run C at differ¬ 
ent times, together with the model-predicted mass function. The horizontal bars on the 
protoplanets are of length 10 th. 
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Fig. 16.— The state of Run C at its endpoint of 10 Myrs, showing eccentricities (top), 
inclinations (middle) and protoplanet masses (bottom) together with their predicted values. 
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Fig. 17.— Snapshots of the protoplanet masses and semimajor axes in Run D at different 
times, together with the model-predicted mass. The horizontal bars on the protoplanets are 
of length 10 Th. 
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Fig. 18.— The state of Run D at its endpoint of 5 Myrs, showing eccentricities (top), 
inclinations (middle) and protoplanet masses (bottom) together with their predicted values. 








Table 1: Simulation parameters 


Quantity 

Run A 

Run B 

Run C 

Run D 

initial solids surface density (g/cm^) 

'“(■riu) 

(rtn) 

(liu) 

225 (riu) 

gas volume density (10“®g/cm^) 


/ \- 13/4 

ii-Hriu) 

/ \ - 13/4 

ii-Hriu) 

10-5 (riu) 

initial radial extent 

1.5 - 5 

5 - 15 

5 - 15 

5 - 15 

initial protoplanet mass (M 0 ) 

0.1 

0.2 

0.2 

0.5 

initial protoplanet number 

22 

16 

16 

12 

planetesimal mass 

0.01 

0.02 

0.02 

0.05 

initial planetesimal number 

4886 

3233 

3233 

3446 

planetesimal gas drag radius 

10 km 

10 km 

1 km 

10 km 





